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ABSIRACT 


The nature of generated machine error in finite digital calculations 
is discussed. The arithmetic of range numbers is developed, and 
examples are given demonstrating the use of range arithmetic as a 
tool for automatic error analysis. A computer program is developed, 
utilizing the TYPE OTHER feature of FORTRAN-63 in conjunction with 
the CDC-1604 digital computer, which enables the user to perform 
automatic error analysis during computation, and a number of programs 


are presented using this feature. 
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l. Introduction. 

The utilization of high-speed digital computers to perform lengthy 
and complex mathematical computations is widespread. Unfortunately, 
the user often assumes that the computer output will necessarily be 
accurate if his input data is accurate and his program logic is sound. 
This need not be the case, and such an assumption may well lead the 
user to incorrect or highly inaccurate results. The difficulty arises 
because of the very nature of digital computation. 

A digital computer has a finite memory, and cannot carry within it 
infinite precision numbers, or perform infinitely accurate calculations. 
Each number stored in a digital computer must, of necessity, contain 
a finite number of digits. As a result, during each step of computation, 
the computer automatically truncates or rounds off numbers to conform 
with its memory storage requirements. Each arithmetic operation 
introduces an error in precision due to such truncation, which may 
accumulate with other such errors to produce sizeable inaccuracies in 
results. It is also possible that the errors may cancel and give a fairly 
accurate result; but, in either case, the user will have no knowledge of 
whether or not any error had accumulated in computation. 

It would be desirable to eliminate precision errors entirely, but 
hardware limitations make it impossible to do so. It is, however, 
possible to estimate and control the errors arising in finite digital 


computations, by appropriate choice of an algorithm. One method is 
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called range arithmetic, or interval arithmetic [1, 2]. This method 
enables the generated error to be analysed automatically within the 
computer, during computation. 

The technique of range arithmetic replaces any real number of infi- 
nite precision by two numbers of finite precision, one of which is a 
lowerbound of the real number, and the other is a corresponding upper- 
bound of the real number. We can then be certain that the real number 
lies in the range of the two finite precision numbers, 

In this thesis, we will discuss the arithmetic of range numbers, the 
application of range arithmetic to digital computation, and present a 
computer program designed for use on the CDC-1604 computer in con- 
junction with FORTRAN-63, which will provide a means of automatic 


error analysis in digital computation. 
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2. The Arithmetik Range amibpers. 

In this section, we will discuss the properties of range numbers, 
and the arithmetic operations pertaining to them. 

Let x be any real number. Associated with this real number is a 
range number X, which has the following properties: 

(1) X isa closed bounded interval; 

(2) X contains x. 

Then define the range number X, associated with the real number x, 
to be the closed bounded interval 

X = XL, xU], 


where xL and xU are defined by the transformations L and U as follows: 
СЭ? 


κε κ xU 
We may then write 
аа X = [xL, xU] 
The transformations L and U are determined by the precision allowed 
in computation. If, for example, we would like to find the range number 
associated with e = 2.71828..... , and we are allowed only two sig- 
nificant figures, then 
Х.-аһдайайабыы R 
On the other hand, if allowed five significant figures, then 


X - [2. 7182, 2. 7183]. 
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For integers (known with exact precision), the range number will have 
the same upper and lower bound, equal to the integer itself. Then, if 
х” Kama}, Sue 

We can see, then, that the range number X will represent the 
smallest closed interval containing x, for a given set of transformations 
L and Un 

The arithmetic operations associated with range numbers can be 
defined in the following manner. Let (op) represent an arithmetic 
operation, and consider two range numbers 


Moai | LS 


and 
Y =[yL, yU]. 
Then, 
ZE X (opp σι 20], 
where 


zL min [xL (op) yL, xL (op) yU, xU(op)yL, xUÚ(op)yU] 


zU = max[xL(op)yL, xL(op)yU, xU(op)yL, xU(op)yU] . 
It can be seen that the real number z - x(op)y will be in the interval Z, 


and that Z represents the smallest interval that contains z. 


Now let 
(+) - denote range addition 
(-) - denote range subtraction 
(*) - denote range multiplication 
(/) - denote range division 
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From the definition above, it can be seen that 


(1) Z X (+) Y = [x ty, zU + ZU) 
(2) -X = [-xU, -xL] 
(3) Z = X(-) Y =-X( ) (` YV) ` [xL -yUÚ?! xU -yE] 
(4) X - Y ifand only if 
xL - yL, and xU - yU 
(5) Z = X(*)Y, where 
zL —= men Γκ Ἑπγο хї ο, ον, xU * JU] 
zZU-- max[xL*yL;, xL &yU, xXU*yL, xU*yU ] 
The values of zL and zU can also be obtained using a table of sign 
discriminations. 


Using the table on page 14, we need only calculate one product for 


each end point, except in sign discrimination 5. 


(6) Z = X(/) Y, where 
ΣΙ, = min [xL /yL; ο το xU το xU / yU ] 
zU = тах [Блуа SEINE απ το] 


Note, that in the case in which Y brackets zero, i.e., Y -[-1, 3], 
range division is undefined. This is equivalent to forbidden division by 


zero in real arithmetic. 
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We now will give some illustrative examples of computations using 


range arithmetic. 


La, b J (1) [ecu [aler bue] 


ра, Баб) Ге с у= [car cb] for c > 0 
[a, b] (*) [c, с] = (сь, са] for c < 0 
Pê) a, b] = [a-b, b-a] 

ΕΤΕ E MEA US fora>0 


ΜΠ. πω» چ‎ 
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3. Range Arithmetic and Automatic Error Analysis. 

In finite digital computation, errors are generated in the process 
of round-off, truncation, and normalization. One has no method of 
determining the seriousness of such errors unless some sort of error 
analysis is performed. A method is necessary that can evaluate the 
error generated at each step in the computation, and also keep track of 
the accumulated error, The techniques of range arithmetic are suited 
to this purpose, If, in addition to doing real arithmetic, we parallel the 
computation with range arithmetic, the range numbers calculated will 
keep track of the error, and the final result will be the smallest closed 
interval containing the ''true"', infinite precision result. 

Error analysis via range arithmetic not only furnishes information 
as to the size of machine-generated errors, but it also can be utilized 
to evaluate the accuracy of alternative computational schemes. Indeed, 
the error generated in one type of computation may be quite different 
from that obtained using another method of computation. A few examples 


will illustrate this point. 


Example 1. 


Consider the problem of finding the roots of the equation 


2 
x + 100,000x + 100 = O . 
We will concern ourselves with finding the larger root of this 


equation, and we will assume that the computer can carry six significant 
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figures, We will first calculate the result using ordinary arithmetic 
and the quadratic formula, and then use range arithmetic in conjunction 
with the quadratic formula, The notation will be that which is equivalent 
to normalized floating-point operations on a computer, For example, 
the number 100, 000 will be denoted by .100000E6, where E6 denotes the 
power of ten multiplying . 100000. 


The equation 


E + bx + c = 0 


has, as its larger root, 


x = (-b + e acre Za 
We then calculate our solution: 


- . 100000E6 + .100000E11 - .400000E3 
„El 


_ -.100000E6 + А/.100000Е11 - .000000Е11 
Е .2EI 


-.100000Е6 + . 100000Е6 
.2El 


. 000000E00 


or 


Now, in order to do the same problem in range arithmetic, we assume 
some uncertainty in the coefficients of our original equation. We then 
replace the original coefficients by range numbers reflecting the un- 

certainty. For example, we will replace 100,000 by [99999.9, 100,001]. 
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The equation to be solved now is 
© 


4 
Г.999999КО, .100001E1 ] x^ + [.999999E5, .100001E6] x 


Now, 


ac 


4ac 


2 
b - 4ac 


N 54 - 4ас = 
-b +y ο΄ - 4ac 


and 


x = 


* [.999999E2, .100001E3] = 0. 


ISSSTE, 


ЕТО, 


РО. 


Г.999998Е2, 
[. 399999ЕЗ, 


' ӨЕ 1 0,.. 


[.999997E10,. 
[.999998Е5, 


[.999998E5, 


[ -.120000E4, 


[-. 120000 E4, 
(. 399999 E3, 


(-. 299994Е1, 


.100001E6]° 


‚100003E11] 


.100001Е1] х (.999999Е2, .100001Е3 1 


. 100003E3 ] 


. 400012E3 ] 


100003E11] - [. 399999E3, .400012E3 } 


100003Е111 


.100002Е61 


.100002Е61 - (.999999Е5, .100001Е6 1 


. 210000E4 ] 


. 210000Е4 ] 
.400012Е3] 


. 525002E1] 


By using range arithmetic, we have kept accurate error bounds and have 


found that our root is in the interval [-2.9, 5.3], or 


Г-.299994Е1, .525002Е11. 


Our result tells us that the true value 


could be anywhere within the calculated interval. However, the interval 
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is wide, and could contain the value zero, which was our answer using 
real arithmetic, We would like to tighten the interval, so that we may 
be more confident of the computer results, The next example will 

illustrate how, with a little knowledge of sources of error generated in 


finite calculations, we can obtain more satisfactory results. 


Example 2. 


Consider the formula for the larger root of a quadratic equation, 


x = (-b + M is - 4ac) / 2a 
and divide numerator and denominator by b so that 

х = (-1+ Jl -s)b/2a 
where 

s = dei 
Expanding the radical in a power series, we obtain 
ATTE OA s 
The equation for x becomes 
х = r ee = жн; —. )b/2a 
Now, for the equation 
2 


x + 100,000x + 100 = 0 


we get 
(.400000E1) (. 100000E1) (. 100000E3) 
С 1000006) 


. 400000 E-7 


И 


s/2 . 200000 E-7 
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a . 200000 E-15 


Biza = (. 100000E6) / (. 200000E1) (. 100000E1) 
= .500000E5 
Then 
x = -(.500000E5) (. 200000E-7 + .200000E-15) 
= -(.500000E5) (. 200000E-7) 
= - ,100000E-2 
ОЕ. 
х = - .001 


We now solve the same problem using range arithmetic, with the 
modified quadratic formula. 
Our equation is 
[.999999E0, .100001E1]x^ 4 [.999999E5, .100001E6]x 
i1 0000909197. 10000113] 2X0... 
In example 2, we found that 


2 
b 


999998210, E О ООЗЕ11] 


4ac о .4 002283| 
therefore, 


.399999E3,  .400012E3 
[.999998E10, . 1000038111 


[.399987E-7, .400013E-7] 


s/2 = [.399987E-7, .400013E-7]/[.200000El, .200000E1] 


[. 199994Е-7, .200007Е-7] 
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2 

s /8 

2 
s/2+s /8 


b/a 


b/2a 


and 


or, 


[.159987E-14, .160009E-14] 
[.199983E-15, .200012E-15] 
[.199994E-7, .200008E-7] 


[.999999E5, .100001E6]/[.999999E0, .100001E1] 


[.999989E5, .100002E6] 


[.499994E5, .500005E5] 


- [.499994E5, .500005E5] [.199994E-7, .200008E-7] 


[ - .100005E-2, -.999958E-3] 


[ - 00100005, -.000999958] . 


Thus, we can see that the method in example 2 generates far less 


error than the straight application of the quadratic formula. Indeed, 


the previous examples illustrate the use of range arithmetic in evaluating 


a given computational technique. 
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4. QRANGE? - Program Description and Usage. 

The TYPE OTHER capability in the FORTRAN-63 compiler for the 
CDC-1604 computer allows the relatively simple implementation of 
range arithmetic for the user. Using the QRANGE7 package, the 
programmer need only declare all floating-point variables to be 
TYPE RANGE?7(3), and follow the usual rules of FORTRAN programming. 
The QRANGE”7 package supplies the necessary subroutines for range 
arithmetic computations, 

Range arithmetic calculations are made carrying a triple of numbers 
(xF, xL, xU), rather than a double as stated previously. The first 
number in the triple is the ordinary floating-point result, as calculated 
using standard floating-point arithmetic. The remaining numbers in the 
triple are the lower and upper range numbers, respectively. Note that 
xL < xF < х7. In any range calculation, tie "real" part of the triple 
is computed last and remains in the accumulator, so that transfers and 
comparisons (such as the IF statement) will follow the same logical 
branches as in ordinary floating-point arithmetic. 

Standard floating-point arithmetic cannot be used on either lower 
or upper range numbers since, in floating-point calculations, round-off 
and normalization are done automatically, and would reduce the accu- 


racy of the range interval. In QRANGE7, the arithmetic performed on 
range numbers is un-normalized, and separate subroutines are provided 


which truncate the lower range numbers and round up by one bit in the 
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least significant position, or truncate, the upper range numbers as 
required. The un-normalized arithmetic is performed by unpacking the 
operands, putting them in un-normalized form, doing the appropriate 
arithmetic, re-normalizing and repacking the result. 

Input/output may be accomplished in one of two ways. The 
programmer may use EQUIVALENCE statements to transfer input / output 
of range numbers as real numbers, where each range number is equiva- 
lent to three real numbers, or he may use the Input / Output for 
Multi- Word TYPE OTHER package [3] and transfer the range numbers 
directly in and out of the computer memory, The following examples 


will illustrate the use of these methods, and the use of QRANGE7. 


29 
TYPE RANGE7(3) A7, B7, X7 
DIMENSION A7(3), AR(9), B?(3) , BR(9) , X7(3, 3), XR(27) 
EQUIVALENCE (A7, AR), (B7, BR), (X7, XR) 
READ 10, ΑΒ, ΒΒ 
10 FORMAT(8F10. 0) 
DO 20 I-1,3 
DO 20 J=1,3 
20 X7(J,I) = A7(I)*B7(J) 
DO 30 151,25, 3 
30 PRINT 40,XR(I), XR(I*1), XR(I42), XR(I), XR(I41), XR(1+2) 
40 FORMAT(//15X,3(E20. 10, 10X)/15X, 3(020, 10X)) 


END 
23 


In this example, we have assumed that the programmer reads in the 


three parts of the range numbers using real floating-point numbers that 


are simultaneously stored as the three parts of the range numbers by 


use of the EQUIVALENCE statement. The next example will illustrate 


the same program, but using the INPUT / OUTPUT for MULTI-WORD 


TYPE OTHER package [ 3]. 


e 


10 


20 


30 


ТУРЕ RANGE?7(3) A7, B7, X7 
DIMENSION A7(3), B7(3), X7(3, 3) 

CALL READ?7 

READ 10,A7,B7 

FORMAT(8F10. 0) 

DO 20 I=1,3 

DO 20 J=1,3 

X7(J,1) = A7(I)*B7(J) 

CALL PRINT7 

PRINT 30, X7(J,I) 

FORMAT(//15X,3(E20. 10, 10X) / 15X, 3(020, 10X)) 


END 
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5. Matrix Inversion Using- Range Arithmetic. 

In order to test QRANGE7, and evaluate the errors generated in 
a particular program, a number of matrix inversions were run. In one 
case, the matrix to be inverted was generated randomly, using the 
RANF(-1) library function of FORTRAN-63. About 15 programs were 
run in this manner. In the second case, the input matrix was specified, 
and the individual elements of the matrix were assumed to be inaccurate. 
Thus, they were inputed as range numbers. Both cases used the 
MATINV2 subroutine provided by the U.S, Naval Postgraduate School 
computer facility, in conjunction with QRANGE7. Two range sub- 
routines had to be used: ABS7, which takes the absolute value ofa 
range number; and Q0Q06700, which complements the range accumulator 
[5]. (A listing of these subroutines appears in Appendix II.) The 


coding for both cases was as follows: 


Case l. Randomly Generated Matrices 
PROGRAM RANMAT 
TYPE RANGE7(3) XMAT, XINV 
DIMENSION XMAT(5, 5), XINV(5, 5), XJUNK(2000) 
READ 2000,N 
DO 90 I=1,N 
90 XJUNK(I) = RANF( - 1) 


DO 10 I-1,5 
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DO 10 J=1,5 
B(I, J) = RANF( - 1)*1000 
10 XMAT(I,J) = B(I, J) 
CALL PRINT7 
DO. 20. J=1, 5 
DO 20 I=1,5 
20 PRINT 2001, ΧΜΑΤΙ(Ι, 7), ΧΜΑΤΙΙ, J) 
CALL MATINV2(XMAT, XINV, 5, 5) 
CALL PRINT7 
DO 30 J=1,5 
DO 30 I=1,5 
30 PRINT 2001, XINV(I, J), XINV(I, J) 
2000 FORMAT(F10. 0) 
2001 КОКМАТ(//15Х,3(Е20.10,10Х)/15Х,3(020,10Х)) 


END 


Case 2, Input Matrix Containing Range Elements 
PROGRAM MMATRIX 
πμ κ Μο πο ο 
DIMENSION A(3, 3), B(3, 3) 
CALL READ? 
DO 10 J=1,3 
DO 10 I=1,3 


10 READ 2000, A(I, J) 
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CALL PRINT? 
DO 20 J-1,3 
DO 20 I=1,3 
20 PRINT 2001, A(I,J),A(I, J) 
CALL MATINV?2(A, B, 3, 3) 
CALL PRINT? 
DO 30 J=1,3 
DO 30 I=1,3 
30 PRINT 2001, Β(Ι, J), Β(Ι, J) 
2000 FORMAT(3F10. 0) 
2001 FORMAT (//15X , 3(E20. 10, 10X) / 15X, 3(020, 10X)) 


END 


The results obtained using random matrices were highly accurate. 
For 5x 5 matrices, seven decimal places of accuracy were obtained, 
and the range of error was small. It might be noted that the conversion 
of the real random matrix to a matrix of range numbers introduced no 
errors, and hence the resulting range intervals on the inverse elements 
were solely a consequence of digital computation. On the other hand, 
the error generated in the second case was quite large, and decimal 
accuracy varied from zero to one, In fact, large errors were generated 
for input coefficients with ranges of approximately 1.0E-2. It was found, 
in both cases, that the error generated increased with increasing 


matrix dimension, This is to be expected, since the number of 
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calculations done in inverting the matrix is increased, For example, a 
12 x 12 random matrix was inverted, and it was found that only four to 
five significant figures could be expected, vice the seven decimal accu- 
racy in the 5 x 5 case, 

One my circumvent the inaccuracies obtained when the input 
elements are range numbers, by the following device [4]. Suppose we 
wish to invert X, whose elements are range numbers. Take Х, to be 
the matrix of mid-points of the elements of X, so that X = x, το 
The elements of E are range numbers of the form[-e, e]. Obtain 


Xo using range arithmetic, and let 


4 z 2 = 2 
Е’ = m lg ES 5; ОЛ Л) : : 


then the series will 


If the elements of E are small compared to Χρ ; 


converge, and one can write 
/ 
X EX MP I. 


The results of a typical run for an inversion of a randomly generated 


5x 5 matrix are shown on the following page. 
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Element Real Part Lower Upper 






(127! 
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4. 6396189281Е01 
8. 6916184500Е02 
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We will present the real part, and the least significant digits of the 


lower and upper range numbers of X-INVERSE. The exponents will not 


be repeated. 


Element Real Part Lower Upper 


TM -9.1422040 728E-04 
(1, 2) 7.419425 11128204 
ΤΕ 
(1, 4) 
(1, 5) 
EN 
(2,25 
(2, 3) 
(2, 4) 
(2, 5) 
(8; 1) 
(8--2) 
(3, 3) 
(3, 4) 
(3, 5) 
(4, 1) 
(4, 2) 
(4, 3 
(4, 4) 
(4, 5) 
(5, 1) 
(5, 2) 
5. 3) 


5 607 
7249 
864 
4688 
517 


ο 
АХ [oo 
~] 
іл i+ 


9.3418389 435Е-05 
-5.779019 51908-05 
6. 2948717 498Е-04 
-1.122627 041Е-03 


e 
~J 
© 


5709 


A 
~J 
-ᾱ 


1. 6441862 954Е-03 944 965 
Або, — 1625103 | 
Ι. 4714936  870E-03 
5. 5237176 832E-05 
1.0620289 855E-03 
1.1383145 460E-03 
2. 3018511 807E-03 
-3. 6993050 965E-03 973 956 
-2. 2307817 966Е-03 962 
6. 4802494 170Е-04 040 338 
-3. 0018154 9738-03 
-5.102677 3200E-05 
4. 0549999 8205-03 825 
6. 5487719 463E-04 477 
1. 3248530 755E-03 742 766 
-1.1561055 861E-04 
3. 6280024 837Е-04 
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6. Conclusions. 

Range arithmetic is a powerful tool, not only for error estimation 
in computations of a scientific or engineering nature, but for the evalu- 
ation and comparison of alternative numerical algorithms as well. The 
QRANGE7 package provides the user with a simple means for such 
analysis. The drawback in the use of QRANGE7, at this point, is that 
computing time is increased, There are many points where the program 
could be increased in efficiency, notably in the Q1Q04770 and Q1Q05770 
subroutines for range multiplication and range division. The author 
feels, however, that even with the increase in computing time, the pay- 
off in the use of range arithmetic is so great that it more than outweighs 
the disadvantage. 

A few function subprograms have been provided for use in conjunc- 
tion with QRANGE7. They are used exactly as one would use the library 
subroutines in FORTRAN, except a "7" replaces the "F" at the end of 
the function name. One must remember to declare these functions 
TYPE RANGE?7(3) before execution. Much work remains to be done in 
this area, and it is hoped that, eventually, a complete library of func- 
tions will be available for use with QRANGE7. 

It might be noted, at this point, that there are errors generated in 
the conversion of decimal-to-binary and binary-to-decimal numbers 
during input/output, and these have not been taken into account. The 


author feels that these errors are negligible compared to the machine 
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4 


generated error occurring during computation. 

The results obtained using QRANGE”7 indicated that the technique of 
range arithmetic can, indeed, be used to provide error information at 
each stage of computation. Unfortunately, the author did not have the 
time to try it out on some of the more common numerical algorithms, 
such as the Runge - Kutta method for solving differential equations, or 
the Newton - Raphson method for determining roots of um order 
polynomials. It is hoped that this will be done in the future, 

In conclusion, it is felt that range arithmetic is of such value that 
it warrants inclusion in future modifications of FORTRAN or other 


algebraic compilers as a standard TYPE. 


m 
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APPENDIX I 


QRANGE? PROGRAM LISTING [5] 


Entry Points: 
Q1Q00770, Q1Q01770, Q1Q02770, Q1Q03770, 
Q1Q04770, Q1Q05770, Q1Q10770, Q1Q00710, 
(31001710, О1002710, О1003710, О1004710, 
(9057105 OBNGLO7L0, Q1Q10170, (01000700, 
Q1Q01700, Q1Q02700, Q1Q03700, Q1Q04700, 
Q1Q05700, Q1Q10700, Q1Q10070, UNPACK7, 
RNGAD7, RNGMU7, RNGDI7, QUIRKL7, QUIRKU7, 


REPACKT, BRROR?777 
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